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We study the phase transition from a nematic phase to a high-density disordered phase in systems 
of long rigid rods of length k on the square and triangular lattices. We use an efficient Monte 
Carlo scheme that partly overcomes the problem of very large relaxation times of nearly jammed 
configurations. The existence of a continuous transition is observed on both lattices for k = 7. 
On the square lattice, our best estimates of the critical exponents differ from those of the Ising 
model, but we cannot rule out a crossover to Ising universality class at larger length scales. On the 
triangular lattice, the critical exponents are consistent with those of the two dimensional three-state 
Potts universality class. We study the correlations in the high-density non-nematic phase, and find 
evidence of a very large correlation length > 2000. 

PACS numbers: 64.60.Cn,64.70.mf,05.50.+q 



I. INTRODUCTION 

The study of the ordering transition in a system of 
long rigid rods in solution with only excluded volume 
interaction has a long history, starting from Onsager's 
realization that they would undergo a transition to an 
orientationally ordered state at large densities In 
two dimensional continuum space, when the rods may 
orient in any direction, the continuous rotational sym- 
metry remains unbroken at any density. However, the 
system undergoes a Kosterlitz-Thouless type transition 
from a low-density phase with exponential decay of orien- 
tational correlations to a high-density phase with a power 
law decay 

In this paper, we study the problem when the underly- 
ing space is discrete: the square or the triangular lattice. 
Straight rods occupying k consecutive sites along any 
one lattice direction will be called fc-mers. For dimers 
(k = 2), it has been shown rigorously that the system 
remains in the isotropic phase at all packing densities 
A system of dimers with additional interactions dis- 
allowin g n earest neighbor occupation can have ordered 
phases For k > 2, the existence of a phase transi- 
tion remained unsettled for a long time [ll|. Ghosh and 
Dhar recently argued that fc-mers on a square lattice, 
for k > k m i n , would undergo two phase transitions, and 
the nematic phase would exist for only an intermediate 
range of densities p* < p < p% [l2| . Similar behavior 
is expected in higher dimensions. In two dimensions, 
numerical studies have shown that k m i n = 7 [l2l ]. The 



'Electronic address: joyjitOimsc.res.m 
tElectronic address: rrajcsh@imsc.rcs 



i.m| 



t Electronic address: ddhar@thcory.tifr. res. in 
§ Electronic address: jstilck@if.uff.br 



r 



existence of the nematic phase, and hence the hrst transi- 
tion from the isotropic to nematic phase, has been proved 
rigorously [l3||. This transition has also been studied in 

detail through Monte Carlo simulations nnni. on the 

square lattice, the transition is numerically found to be 
in the Ising 11411 . or equivalently in the liquid-gas univer- 
sality class |l8j |. and on the triangular lattices, it is in the 
q = 3 Potts model universality class [3, EH ■ 

In this paper, we ask whether the high-density disor- 
dered phase is a reentrant low-density disordered phase, 
or a new qualitatively distinct phase. To distinguish be- 
tween these two phases without nematic order, we will 
refer to the first as low-density disordered (LDD) phase 
and the second as high-density disordered (HDD) phase 
in the remainder of the paper. 

The second transition at p\ from the nematic to the 
HDD phase has not been studied much so far. Numerical 
studies arc difficult because of the large relaxation times 
of the nearly jammed configurations at high densities. 
Conventional Monte Carlo algorithms using deposition- 
evaporation moves involving only addition or removal of 
single rods at a time are quite inefficient at large den- 
sities. With additional diffusion and rotation moves, it 
is possible to equilibrate the system [13, EH , but the al- 
gorithm is still not efficient enough to make quantitative 
studies of the transition, or the nature of the HDD phase. 
In Ref. [l2j], a variational estimate of the entropy of the 
nematic and the HDD phases suggests that 1 — p* 2 should 
vary as 1/fc 2 for large k. Linares et. al. estimated that 
0.87 < p* 2 < 0.93 for k = 7, and proposed an approx- 
imate functional form for the entropy as a function of 
the density |17| . However, not much is known about the 
nature of transition. Recently, we showed that a Bethc- 
like approximation becomes exact on a random locally 
tree like layered lattice, and for the 4-coordinated lat- 
tice, k m i n = 4, but on this lattice, the second transition 
is absent 12011 . 
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In this paper, by implementing an efficient Monte 
Carlo algorithm, we show that at high densities the ori- 
cntational order is absent. We investigate the nature of 
this HDD phase. We find evidence of a power-law decay 
of orientational correlations between rods for distances 
r < £* ~ 2000, where £* is a characteristic length scale 
of the system. Interestingly, there are indicators that the 
nature of correlations changes for distances r > £*, but 
we can not rule out a power-law decay of correlations 
even for r 3> £*, as this regime is beyond our present 
computational resources. 

Regarding the critical behavior near the phase tran- 
sition on the square and triangular lattices, for k = 7, 
our results show that the transition is continuous and 
occurs for p* 2 = 0.917 ± .015 on the square lattice, and 
for p\ = 0.905 ± .010 on the triangular lattice. On the 
square lattice, our best estimates of the effective crit- 
ical exponents differ from the Ising universality class, 
with exponents v = 0.90 ± 0.05, fi/v = 0.22 ± 0.07, 
j/u = 1.56±0.07 and ajv = 0.22±0.07. However, it ap- 
pears that these are only effective exponents, and we can- 
not rule out a crossover to the Ising universality class at 
larger length scales. On the triangular lattice, the second 
transition is indistinguishable from the three-state Potts 
model universality class with critical exponents v = 5/6, 
P/v = 2/15, 7/zy = 26/15 and a/v = 2/5. 

The plan of the paper is as follows. In Sec. [Til wc de- 
fine the model precisely, and describe the Monte Carlo 
algorithm used. In Sec. IIII1 we use this algorithm to 
show that at high activities, the nematic phase is un- 
stable to creation of bubbles of HDD phases, and that 
the decay of the nematic order parameter to zero is well- 
described quantitatively by the classical nucleation the- 
ory of Kolmogorov-Johnson-Mehl-Avrami. In Sec. IIV1 
wc study different properties of the HDD phase: the two 
point correlations, cluster size distributions, susceptibil- 
ity, size distribution of structures that we call 'stacks', 
and the formation of bound states of vacancies. The 
critical behavior near the second transition from the ne- 
matic phase to the HDD phase is studied in Sec. [V] for 
both the square and triangular lattices, by determining 
the numerical values of the critical exponents. Section IVT1 
summarizes the main results of the paper, and discuss 
some possible extensions. 



II. MODEL AND THE MONTE CARLO 
ALGORITHM 

For simplicity, we first define the model on the square 
lattice. Generalization to triangular lattice is straight- 
forward. Consider a square lattice of size L x L with 
periodic boundary conditions. A fc-mer, a straight rod 
occupying k consecutive lattice sites in any one lattice 
direction, can be either horizontal (x-mer) or vertical (y- 
mer). A lattice site can have at most one fc-mer passing 
through it. An activity z = is associated with each 
fc-mer, where /x is the chemical potential. 



The Monte-Carlo algorithm we use is defined as follows 
(this was reported earlier in a conference (2~l|): given a 
valid configuration, first, all x-mers are removed. Each 
row now consists of sets of contiguous empty sites, sepa- 
rated from each other by sites occupied by y-mcrs. The 
re-occupation of these empty intervals using x-mers can 
be done independently in each row, and the problem 
reduces to that of occupying an interval of some given 
length £ of a one dimensional lattice with k-mers with cor- 
rect probabilities corresponding to the equilibrium grand 
canonical ensemble. 

Let the partition function on a one dimensional lattice 
of £ sites with open boundary conditions be denoted by 
f2 (£). The probability that the left most site is occupied 
by the left most site of a x-mer is z£l {£ — k)/tl (£). If 
not occupied, we consider the neighbor to the right and 
reduce the number of lattice sites by one. If occupied, we 
move to the (A: + l)* /l neighbor and reduce the number of 
lattice sites by k. 

The partition functions £l (£) obeys the simple recur- 
sion relation Cl (£) = zfl (£—k) + £l (£—l), for £ > k, and 
£l (£) = 1 for £ = 0, 1, . . . , k — 1. The recursion relation 
is easily solved by £l (£) <x \ l . With periodic boundary 
conditions, the recursion relations have to be modified. 
Let Q p (£) be the partition function of a one dimensional 
lattice of length £ with periodic boundary conditions. It 
is easy to see that Sl p (£) = kzQ D (£ — k) + Q a (£ — 1). Wc 
use a list of stored values of the relevant probabilities to 
reduce the computation time. 

Following the evaporation of and re-occupation by x- 
mers, all y-mers are evaporated and all columns reoc- 
cupicd with y-mcrs. A Monte Carlo move corresponds 
to one set of evaporation and re-occupation of both x- 
mers and y-mers. It is straightforward to see that the 
algorithm is ergodic, and satisfies the detailed balance 
condition. 



III. METASTABILITY OF THE NEMATIC 
PHASE FOR LARGE ACTIVITIES 

We first verify that, for large activities, the nematic 
phase is unstable to the growth of the HDD phase. In 
Fig.QJa)-(c), we show snapshots of the system of rods of 
length k = 7 in equilibrium, on a square lattice at low, 
intermediate and high densities. For the high-density 
snapshot, the initial configuration had full nematic order, 
but the system relaxed to a disordered phase. A similar 
disordered phase is also seen for the triangular lattice at 
high densities (see Fig. [2j . 

In Fig. [3l we show the temporal evolution of the or- 
der parameter Q, defined by Q — (n^ — n ll )/(n/ l + n v ), 
where and n v are the number of x-mers and y-mers 
respectively. For all values of /i, the initial configuration 
had full nematic order. For /i = 3.89, at large times, 
the system relaxes to an equilibrium state with a finite 
nematic order. However, for larger /.t = 7.60, the ne- 
matic order decreases with time to zero. Interestingly, 
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FIG. 1: Typical configurations of the system in equilibrium at densities (a) p m 0.66 (p = 0.41) (b) p ~ 0.89 (p = 4.82), and 
(c) p ~ 0.96 (p = 7.60) on a square lattice. Here, k = 7 and L = 98. 
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FIG. 2: A typical configuration of the system in equilibrium 
at density p m 0.96 (/i = 7.60) on a triangular lattice. Here, 
k = 7 and L = 98. 



we find that the average lifetime of the metastable state 
decreases with increasing system size, and saturates to a 
L independent value for L > 200 (sec Fig. [3]). 

Naively, faster relaxation for larger systems may ap- 
pear unexpected, but is easily explained using the well- 
known nucleation theory of Kolmogorov-Johnson-Mchl- 
Avrami (2^, HH . We assume that critical droplets of the 
stable phase are created with a small uniform rate e per 
unit time per unit area, and once formed, the droplet 
radius grows at a constant rate v. Then, the probability 
that any randomly chosen site is still not invaded by the 
stable phase is given by exp[-e f* dt'V(t% where V(t') 
is the area of the region such that a nucleation event 
within this area will reach the origin before time t' . The 
area V(t') is given by V(t') = Trv 2 t' 2 when the droplet is 
smaller than the size of the lattice. For time t' greater 
than this characteristic time t* , we have V(t') = L 2 . If 



FIG. 3: Decay of the order parameter Q for the square lat- 
tice as a function of time (Monte Carlo steps), starting from 
a fully ordered state for two different values of p: p, — 3.89 
(p » 0.867), and p = 7.60 (p w 0.957). The best fit of the data 
to Eq. {TJ with additional subleading terms is also shown. In- 
set: Data for different chemical potentials, all corresponding 
to HDD phase for L = 154 and k = 7. The densities corre- 
sponding to these values of p are approximately 0.957, 0.948, 
0.941. 



the droplet does not grow equally fast in all directions, 
we take suitably defined average over directions to define 
v 2 . Thus, we obtain 



Q(t) 



cxp 
exp 



-irevH* 2 | t 



, for t < £*, 

2t 
~3 



for t>t*. (1) 



We see that with this choice, both Q{t) and its deriva- 
tive are continuous at t = t*. Since V(t') should tend to 
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L 2 for large t', we get the crossover scale t* given by 



t* = 



L 



(2) 



The crossover lattice size L* beyond which the average 
lifetime of the mctastable state becomes independent of 
L can then be estimated from the above to be 



L* 



/3y^ 1/3 



(3) 



Fitting the numerical data in Fig. [3] to Eq. (JXJ) we ob- 
tain e = (2.1 ± 0.2) x 1(T 10 and v = (5.5 ± 0.7) x 10~ 5 
for \i = 7.60. From Eq. ([3]), we then obtain the crossover 
scale L* ~ 110, of the same order as the numerically ob- 
served value of L* ~ 200. The difference is presumably 
due to simplifications made in the theory, e.g., neglect- 
ing the dependence of the mean velocity of growth on the 
direction of growth, or the curvature of the interface, etc. 

We can also estimate v directly from simulations of a 
system with an initial configuration where half the sam- 
ple is in the nematic phase and the other half is in the 
equilibrium disordered phase at that /x. For /i = 7.60, we 
find that the velocity reached its asymptotic L indepen- 
dent value slowly, and is equal to 10 x 10 -5 for L = 784, 
reasonably close to the velocity obtained from fitting data 
to Eq. ([T]). For decreasing chemical potential fi, we find 
that both the velocity v and nuclcation rate e increase. 



IV. NATURE OF THE HIGH-DENSITY 
DISORDERED PHASE 

Based on the behavior seen in other hard-core systems, 
we expect the behavior of different correlations in the 
HDD phase to be qualitatively different from the known 
exponential decay of correlations in the LDD phase. For 
a system of oriented hard squares, or hard rectangles, 
in continuum space, one expects the high-density phase 
to be qualitatively different from the low-density phase: 
for example, oriented hard squares are expected to show 
columnar order, or quasi-long-range bond-oricntational 
order. This suggests that HDD phase for the fc-mer prob- 
lem is also not the same as the LDD phase. In this sec- 
tion, we test this possibility by studying the susceptibil- 
ity Xi the order parameter correlation function Cqq(%, j), 
the cluster size distribution F(s), and the size distribu- 
tion of structures that we call stacks. We also examine 
the formation of bound states of vacancies. 

The susceptibility is defined as x = L 2 ((rih — 
n v) 2 ) I (nh + n v ) , where nu and n v are the number of 
x-mcrs and y-mcrs. Figure 0] shows the variation of x 
with L, for three different values of fi in the HDD phase. 
There is no anomalous dependence on L, hence, if the 
correlations are a power law, then the decay exponent is 
larger than 2. This implies that the order parameter Q 
should scale as which is confirmed in the inset of 
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FIG. 4: Susceptibility % f° r the square lattice as a function 
of L for three values of fi, all in the HDD phase. There is no 
anomalous dependence on L. Inset: The scaled probability 
distribution for the order parameter P(Q) for different L's 
collapse when plotted against QL. 



Fig. SI where the scaled probability distributions for dif- 
ferent L's collapse onto one curve when plotted against 
QL. 

The order parameter correlation function Cqq(i,j) is 
defined as follows. Given a configuration, wc assign to 
each site a variable Sij, where Sij = 1 if is 

occupied by an x-mer, Si^ = — 1 if is occupied by 
an y-mer, and Sij = if (i, j) is empty. Then, 



Cqq(*7 j) = (SoflSij). 



(4) 



Figure [5] shows the variation of Cqq (r) with separation 
r along the x- and y- axes, for different chemical poten- 
tials and systems sizes. In the HDD phase, the correla- 
tion function has an oscillatory dependence on distance 
with period k, and for r ^ k, appears to decrease as a 
power law r _?) , with r\ > 2. Given the limited range of r 
available 7 « r < L/2, it is difficult to get an accurate 
estimate of the exponent rj. 

The long-range correlations in the HDD phase are bet- 
ter studied by looking at the large-scale properties of con- 
nected clusters of parallel rods. For instance, it is known 
that the exponent characterizing the decay of cluster size 
distribution of critical Fortuin-Kasteleyn clusters 24j in 
the g-state Potts model [1^, [26| has a non-trivial depen- 
dence on q. We denote all sites occupied by x-mers by 
1 and the rest by zero. For our problem, we define a 
cluster as a set of l's connected through nearest neigh- 
bor bonds. Let F (s) be the probability that a randomly 
chosen 1 belongs to a cluster of s sites. Clearly, F(s) 
is zero, unless s is a multiple of k. Let the cumulative 
distribution function be F cum (s) = X) s '=i ^( s ')- 

In Fig. [6l we plot F cum (s) in the HDD phase for dif- 
ferent system sizes. We find that for intermediate range 
of s, for 10 3 < s < 10 6 , F cum (s) ~ As 1 "", with r < 1. 
For fi = 7.60, we estimate the numerical values to be 
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FIG. 5: Order parameter correlations CQQ(r) for the square 
lattice as a function of r, measured along the x- and y- axes, 
for three different values of fi, all of them larger than /i c ~ 
5.570. The system size is L — 252. Inset: The dependence of 
CQQ(r) on L is shown for fj, — 7.60. The solid lines are power 
laws r~ 2 ' 5 , intended only as guides to the eye. 
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FIG. 6: Fcum(s), the probability that a randomly chosen 1 
(a site occupied by a x-mer) belongs to a connected cluster 
of size < s, in the HDD phase (fj, = 7.60) for different system 
sizes. The data are for the square lattice. 



A = 0.037 and t = 0.762. For small system sizes (up to 
L = 1568), F cum (s) has a system-size dependent cutoff. 
The L-independent cutoff s* is determined by the con- 
dition As* 1 " 7 ' w 1, giving s* w 1.04 x 10 6 . The density 
of l's being roughly 0.48, we expect to observe s* only 
when L exceeds a characteristic length scale £* ~ 2000. 
This is indeed seen in Fig. [6l 

In the HDD phase, F cum (s) depends weakly on fj, (see 
Fig. [71) ■ The power law exponent r is estimated to be 
0.778 (fj, = 6.50), 0.767 (fj, = 6.91) and 0.762 (fi = 7.60). 
It appears that r increases slowly with increasing [i^ while 
s* decreases with increasing [i. 

One qualitative feature of the HDD phase is the ap- 
pearance of large groups of parallel rods, worm-like in 



FIG. 7: Fcum(s), the probability that a randomly chosen 1 
(a site occupied by a x-mer) belongs to a connected cluster 
of size < s, for different [i, all corresponding to the HDD 
phase. The curves appear to have weakly density dependent 
power-law exponents. 




FIG. 8: Some examples of the different types of stacks, shown 
here as rods joined by wiggly lines, for (a) square lattice and 
(b) triangular lattice. The snapshots are for fi = 7.60, cor- 
responding to the HDD phase. Rods of different orientations 
are shown in different colors for easy visualization. 



appearance, nearly aligned in the transverse direction. 
This is clearly seen in Fig. (He). Wc call these groups 
stacks. To be precise, we define a stack as follows: two 
neighboring parallel fc-mers are said to belong to the same 
stack if the number of nearest-neighbor bonds between 
them is greater than k/2. A stack is the maximal cluster 
of rods that can be so constructed. By this definition, a 
stack has a linear structure without branching, with some 
transverse fluctuations allowed. Examples of stacks on 
square and triangular lattices are shown in Fig. [51 Any 
given configuration of rods is uniquely broken up into a 
collection of disjoint stacks. 

There are a noticeable number of large stacks in the 
HDD phase. We measured the stack size distribution 
D(s), the number of stacks of size s per site of the lattice, 
in all the three phases and at the transition points (see 
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FIG. 9: Stack distribution in the LDD phase (jj, = 0.200), 
intermediate density nematic phase (p, = 3.476), HDD phase 
(/i = 7.600), and at two critical points (p — 1.3863, 5.570) are 
shown. Data are for L — 280, k = 7, and the square lattice. 
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FIG. 10: The average distances da and dij, on the square 
lattice, between a vacancy on sublattice i and the nearest 
vacancy on sublattice j as a function of density p. The i 7^ j 
data are averaged over all j ^= i. The solid lines show the 
functions K(l - p)~ 1/2 , for K = 1.36 and 1.12. 



Fig. IH]). Interestingly, we found that this distribution 
is nearly exponential in all the three phases, as well as 
at the critical points, and there is no indication of any 
power-law tail in this function. In the HDD phase, the 
mean stack size is approximately 12, for both square and 
triangular lattices, and is only weakly dependent on the 
density. 

It was suggested in [l| that the second phase transi- 
tion may be viewed as a binding-unbinding transition of 
k species of vacancies. For studying such a characteriza- 
tion, we break the square lattice into k sublatticcs. A site 
(x, y) belongs to the i-th sublattice if x + y = i (mod k), 



where 



0, 1, . . . , k — 1. Let dij be the mean value of 



the Euclidean distance between a randomly picked va- 
cant site on the i-th sublattice, and the vacant site near- 
est to it on the j-th sublattice. In a typical configura- 
tion with a low density of vacant sites, we would expect 
the vacancies to form bound states with k vacancies per 
bound state, one vacancy on each sublattice. The HDD 
phase can be described weakly interacting gas of 
such bound states, when the typical distance between 
two bound states is much larger than the mean size of a 
bound state i.e. da d^, for j =/= i. 

In Fig. 1101 we show the variation of and da with 
density p. We see that da and dij with i ^ j both vary 
approximately as (1 — p) -1 / 2 , with da w 1.2dij. The data 
are for L = 168 and k = 7. There is no noticeable depen- 
dence of the data on L. We see that the saturation value 
of dij (i 7^ j) is not reached for the largest densities that 
we simulated. We conclude that size of the bound state, 
if it exists, would be substantially greater than 7 lattice 
units. Near p\ 1 the typical spacing between vacancies 
is much less than the the size of the bound state, and 
the transition can not be understood as a density-driven 
binding-unbinding transition of vacancies. 



V. CRITICAL BEHAVIOR NEAR THE SECOND 
TRANSITION 

We now discuss the critical behavior near the sec- 
ond transition. Several thermodynamic quantities are 
of interest. We define the nematic order parameter m 
as follows. For the square lattice, m = (n/j — n v )/N, 
where i%h and n v are the number of lattice sites occu- 
pied by x-mers and y-mers respectively, and N is the 
total number of lattice sites. For the triangular lattice, 
m = (ni +w«2 +uj 2 n3)/N , where ui is the complex cube- 
root of unity, and 711,712,713 are the number of sites oc- 
cupied by fc-mcrs oriented along the three allowed direc- 
tions. The density p is defined by the fraction of sites 
that are occupied by the fc-mers. The order parameter 
Q, its second moment Xi compressibility n, and Binder 
cumulant U are defined as 



Q 



(H) 
(p) ' 

L 2 (\m\*) 
(P) 2 ' 
L 2 [(P 2 ) - (P) 2 ] 



U 



a(\m\ 



(5a) 

(5b) 
(5c) 
(5d) 



where a = 3 for square lattice and a = 2 for triangu- 
lar lattice. Q is zero in the LDD and HDD phases and 
nonzero in the nematic phase. 

The data used for estimating the critical exponents 
are for k = 7, and for five different system sizes L = 154, 
210, 336, 448, and 952 for the square lattice and L = 210, 
336, 448, and 560 for the triangular lattice. The system is 
equilibrated for 10 7 Monte Carlo steps for each /Lt, follow- 
ing which the data are averaged over 3 x 10 8 Monte Carlo 
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FIG. 11: The temporal variation of the autocorrelation func- 
tion AQQ(t), where Aqq(0) has been normalized to 1. Q rep- 
resents (a) the system order parameter and (b) the local single 
site order parameter, where the site takes values 1,-1 depend- 
ing on whether a x-mei or y-mer passes through it, otherwise 
being zero. The data are for ft = 7.60 and the autocorrelation 
times corresponding to the solid lines are (a) 220000 and (b) 
52000. Inset: Data as above but in the LDD phase (fi — 0.20), 
with (c) corresponding to the system order parameter and (d) 
corresponding to the local single site order parameter. The 
autocorrelation times are (c) 440 and (d) 310. All data are 
for k = 7. 



steps. These times are larger than the largest autocorre- 
lation times that we encounter (see Fig. [TT]) . The largest 
autocorrelation time is for the largest density and is close 
to 2.2 x 10 5 Monte Carlo steps. To estimate errors, the 
measurement is broken up into 10 statistically indepen- 
dent blocks. We check that the long time behavior of the 
system is independent of the initial condition by start- 
ing from two different initial conditions, one in which all 
sites are occupied by x-mers and the other in which one 
half of the lattice contains x-mers and the other half only 
y-mers. 

The quantities in Eq. (|5|) are determined as a function 
of fi using Monte Carlo simulations. The nature of the 
second transition from the ordered nematic phase to the 
HDD phase is determined by the singular behavior of U, 
Q, x, and k near the critical point. Let e = (fi — fi c )/ fj, c , 
where fi c is the critical chemical potential. The singular 
behavior is characterized by the critical exponents v, /?, 
7, and a, defined by Q ~ (— e)' 3 , e < 0, x ~ N~ 7 an( i 
K ~ |e|~ Q , and £* ~ | e | ^ , where £* is the correlation 
length and |e| — > 0. The exponents are obtained by finite 
size scaling of the different quantities near the critical 
point: 

U ~ MeL 1 ^), (6a) 
Q ~ L-^f^eL 1 /"), (6b) 
X * L-y/»f x (eL^), (6c) 
k ~ L a ' v f K {eL 1 ' v ), (6d) 



4.8 5 5.2 5.4 5.6 5.8 6 

M 

FIG. 12: The Binder cumulant U as a function of chemical 
potential fi for different lattice sizes of a square lattice. The 
curves intersect at fi c = 5.570 ± .02. Inset: Data collapse for 
square lattices when U is plotted against eL 1 ^ with v = 0.90 
and e = (fj, — fJ, c )/fx c - 



where /„, f q , f x , and f K are scaling functions. 



A. Square lattice 

We first present results for the square lattice. The 
data for the Binder cumulant U for different system sizes 
intersect at [i c = 5.570 ± .02 (sec Fig. [15]). The density 
at this value of chemical potential is p* 2 = 0.917 ± .015, 
consistent with the variational estimate 0.87 < p* 2 < 0.93 
in Rcf. [ItJ ■ The value of U at the transition lies in 
the range 0.56 to 0.59. This is not very different from 
the value for the Ising transition U c ~ 0.61 [13]. The 
data for different system sizes collapse when scaled as in 
Eq. ([6aD with v = 0.90 ± .05 (see inset of Fig. [12]). To 
compare with the first transition from the LDD phase to 
the nematic phase, p\ = 0.745±0.010, and the exponents 
are consistent with the Ising model exponents, so that 
v= 1 P. 

The data for order parameter, x an d n for different sys- 
tem sizes are shown in Figs. [TjJ] [T3] and Q1)] respectively. 
Q decreases to zero at high densities. Our best estimates 
of effective critical exponents are f3/is = 0.22 ± 0.07 (see 
inset of Fig. [13]). j/u = 1.56 ±0.07 (see inset of Fig. [H]), 
and a/v = 0.22 ± 0.07 (see inset of Fig. [15]). The es- 
timated error bars are our subjective estimates, based 
on the goodness of fit. These differ substantially from 
the values of the exponents of the two dimensional Ising 
model {y = 1, ft = 1/8,7 = 7/4, a = 0). However, as 
discussed in Sec. IIVI it seems like there is a character- 
istic length scale £* ~ 2000 in the HDD phase, and we 
cannot say much about the asymptotic value the critical 
exponents at length scales L ^> £*. 
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FIG. 13: The variation of the order parameter Q with chem- 
ical potential fi for different systems sizes of a square lat- 
tice. Inset: Data collapse for square lattices when scaled 
Q is plotted against eL 1/l/ with v = 0.90, fi/v = 0.22 and 
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FIG. 14: The variation of \, the mean of the square of the 
order parameter, with chemical potential fi for different sys- 
tem sizes of a square lattice. The curves cross at fi c when \ 
is scaled by L _7// ", with j/v = 1.56. Inset: Data collapse 



for square lattices when \L 
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plotted against eL 1/V with 



v = 0.90, and e = (n — \i^)/\i c 



FIG. 15: The variation of compressibility k with chemical 
potential /j, for different system sizes of a square lattice. Inset: 
Data collapse for square lattices when the scaled k is plotted 
against eL 1 ^ with v = 0.90, a/u = 0.22, and e = (fi — f± c )/n c . 
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FIG. 16: Finite size scaling of the triangular lattice data of (a) 
U, (b) Q, (c) x, and (d) k as in Eq. © with u = 5/6, /3/u — 
2/15, y/u = 26/15 and a/u = 2/5. The critical chemical 
potential is [i c — 5.147. 



B. Triangular lattice 

For the triangular lattice, we find that the second tran- 
sition is continuous with /i c = 5.147 ± .05 and p\ = 
.905 ± .005. The data for U, Q, x> an( i K f° r different 
system sizes collapse onto one scaling curve when scaled 
as in Eq. $Q with exponents that are indistinguishable 
from those of the three state Potts model (see Fig. fTB|) 
= 5/6, = 1/9, 7 = 13/9 and a = 1/3). However, 
we find that in the HDD phase, like in the square lattice, 
the Q-Q correlations along the three directions of fc-mer 
orientations decrease as a power law for the lattice sizes 



that we have gone up to. 



VI. SUMMARY AND DISCUSSION 

In this paper, we studied the problem of hard, rigid 
rods on two dimensional square and triangular lattices, 
using an efficient algorithm that is able to overcome jam- 
ming at high densities. We showed the existence of a 
second transition from the ordered nematic phase to a 
disordered phase as the packing density is increased. By 
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studying the order parameter, its second moment, com- 
pressibility and Binder cumulant, we concluded that the 
second transition is continuous on both square and tri- 
angular lattices. We also investigated the nature of cor- 
relations in the HDD phase by measuring distribution of 
connected clusters of parallel rods, as well as the distri- 
bution of stacks. 

An interesting feature of the HDD phase is the appear- 
ance of a large correlation length scale of £* ~ 2000 on the 
square lattice, as inferred from the fact that the cluster 
size distribution seems to follow a power law distribution 
F(s) ~ S ~ T , with r < 1 for s < £* 2 . The amplitude 
of this power-law term is rather small. This is related 
to the fact that for the fc-mer problem, various pertur- 
bation series involve terms like k~ k , which then leads to 
large correlation lengths. The HDD phase has power-law 
correlations at least for lengths up to £*, perhaps for all 
r. 

On the square lattice, our best estimates of the nu- 
merical values of the critical exponents are different from 
those of the Ising universality class. However, because 
the correlation lengths in the HDD phase are large, we 
cannot rule out a crossover to the Ising universality class 
at larger length scales. For the triangular lattice with 
k = 7, the estimated exponents for second transition are 
indistinguishable from those of two dimensional 3-state 



Potts universality class with critical exponents v = 5/6, 
P/v = 2/15, 7/1/ = 26/15 and a/v = 2/5. 

We expect that the nature of the transition will not 
depend on the rod length k as changing k is equivalent 
to change in the length scale. It has been predicted that 
p* = 1 — a/k 2 + . . . for large k, where a is a constant 
[l2| . It would be interesting to verify this conjecture by 
studying the transition for larger k on both triangular 
and square lattices. 

The Monte Carlo algorithm that we have implemented 
is more efficient than Monte Carlo algorithms with local 
single particle moves. In addition to overcoming jamming 
at high packing densities, it is easily parallelized, which 
makes it suitable for studying hard-core systems with 
other particle shapes, and also in higher dimensions. This 
seems to be a promising approach for future studies. 
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